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ABSTRACT 


As  lasers  become  more  pervasive,  the  dangers  posed  by  laser  radiation  to  the  human  retina  have 
increased.  Computer  aided  modeling  of  laser  tissue  interaction  for  energies  deposed  in  the  retina 
allows  for  researchers  to  simulate  parameter  ranges  where  exposure  can  cause  damage.  While 
this  approach  has  been  studied  using  a  variety  of  methods,  one  of  the  most  referenced  models  has 
been  the  Thompson-Gerstman  melanin  granule  model.  Due  to  limited  computing  resources  at 
the  time,  the  original  FORTRAN  version  of  this  model  was  implemented  as  a  stand-alone  serial 
code  and  was  only  able  to  model  single-pulse  exposures.  A  new  C++  version  of  the  Thompson- 
Gerstman  model  has  been  implemented,  which  expands  both  the  functionality  and  portability  of 
the  method.  The  source  type  has  been  expanded  to  include  multi-pulse  as  well  as  single  pulse. 
Additionally,  the  functions  developed  by  Thompson  and  Gerstman  are  now  part  of  a  library  that 
can  be  accessed  through  other  modeling  software. 
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1  INTRODUCTION 


Laser  radiation  in  the  visible  and  near  infrared  (400-1400  nm)  region  of  the  spectrum  poses  a 
unique  risk  to  the  human  eye.  A  collimated  laser  beam  incident  on  the  cornea  will  be  focused  to 
a  small  spot  on  the  retina  with  a  transmission  between  the  cornea  and  retina  of  approximately 
90%  at  550  nm  [1],  As  light  passes  through  the  retina,  the  layers  of  tissue  will  absorb  energy, 
which  could  lead  to  damage.  Damage  to  the  retina  most  likely  occurs  in  the  retinal  pigmented 
epithelium  (RPE)  layer  where  50%  of  600  nm  light  is  absorbed  [2],  RPE  cells  in  the  pigment 
layer,  located  at  the  back  of  the  retina,  may  be  irreparably  damaged  if  laser  radiation  deposition 
is  sufficiently  intense.  Damage  to  cells  in  the  RPE  layer  may  result  in  loss  of  vision  sensing  in 
the  affected  region  of  the  retina. 

Damage  can  be  caused  by  three  distinct  mechanisms:  photo  thermal,  photomechanical,  and 
photochemical.  Photothermal  damage  occurs  when  the  temperature  of  tissue  rises  to  levels 
where  protein  denaturation  or  enzyme  inactivation  occurs.  This  type  of  damage  is  produced  with 
a  long-pulse  duration  (t  >  10  ps)  where  substantial  energy  is  deposited  slowly.  Photomechanical 
damage  occurs  under  short  and  ultrashort  pulses  (10  ps  <  t  <  10  ps)  where  energy  is  deposited 
quickly,  causing  mechanical  stress  or  micro-cavitation  within  the  tissue  and  cells. 

An  accurate  computational  model  of  the  laser-tissue  interaction,  which  accounts  for  energy 
absorption  and  thermal  diffusion,  could  be  used  to  predict  damage  thresholds  based  on  the 
properties  of  the  laser  source.  The  relationship  between  beam  characteristics  such  as 
wavelength,  pulse  duration,  beam  profile,  and  spot  size  could  then  be  compared  to  damage 
thresholds  and  lesion  size,  which  would  give  insight  into  exposure  thresholds  in  regimes  where 
no  laser  damage  studies  have  been  conducted.  This  tool  would  have  the  added  benefit  of 
reducing  the  dependence  on  animal  studies  for  detennining  damage  thresholds. 

The  topic  of  modeling  photothennal  injury  to  the  retina  has  been  evaluated  by  a  wide  range  of 
authors  [3-8].  While  several  possible  simulations  have  emerged  to  model  the  laser-tissue 
interaction  in  the  retina,  a  “granular  model”  has  advantages  over  other  models,  which  try  to 
numerically  solve  the  heat  equation  assuming  homogeneous  retinal  layers  and  homogeneous  heat 
absorption.  At  visible  and  near-infrared  wavelengths  most  optical  absorption  in  the  RPE  is 
within  the  intracellular  melanin-impregnated  granules  known  as  melanosomes.  For  pulses 
shorter  than  0.1  ms,  the  assumption  of  homogeneous  temperature  within  layers  breaks  down,  as 
heat  near  the  energy  absorbing  melanosomes  is  much  greater  than  in  surrounding  tissue. 

A  number  of  researchers  have  developed  granular  models  [9-13]  in  order  to  address  the  issue  of 
non-homogeneous  heating  of  the  RPE.  The  most  advanced  of  the  granular  models,  and  the  one 
most  often  referenced  by  other  researchers  since  its  development,  is  the  Thompson-Gerstman 
melanin  granule  model  [14].  Instead  of  solving  for  temperature  across  the  domain  through  finite 
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element  or  finite  difference  meshing,  the  Thompson-Gerstman  model  instead  solves  the  analytic 
heat  profile  of  a  single  melanosome,  then  calculates  the  temperature  at  any  given  point  in  the 
medium  through  superposition  of  thermal  energy.  This  report  documents  the  creation  of  a  C++ 
implementation  of  the  Thompson-Gerstman  model  by  Austin  Kane,  which  replaces,  and 
improves  on,  the  original  FORTRAN  implementation  created  in  1996  by  Dr.  Randy  Thompson. 
The  C++  implementation  should  allow  for  greater  utility  and  easier  integration  into  other  laser- 
tissue  damage  models. 


2  MODELING  LASER-TISSUE  INTERACTION 
2.1  Thompson-Gerstman  Model 

An  analytic  expression  for  heated  spheres  in  a  homogeneous  infinite  medium,  originally  derived 
by  Goldenburg  in  the  1950’s  [9,10],  forms  the  basis  of  the  Thompson-Gerstman  granular  model. 
A  granular  model  can  also  be  modified  to  include  photomechanical  damage  from  bubble 
formation  near  the  granules,  which  has  been  presented  in  other  models  [11].  However,  the 
Thompson-Gerstman  model  considers  only  photo  thermal  effects.  While  many  granule  models 
assume  melanosomes  of  zero  diameter  to  reduce  the  complexity  of  the  analytical  expression,  this 
assumption  breaks  down  for  ultra-fast  pulses  where  calculating  the  temperature  rise  near 
granules  becomes  inaccurate  [12].  The  Thompson-Gerstman  model  has  melanosomes  of  finite 
size,  where  the  analytical  solution  of  the  heat  equation  has  been  derived  at  the  granule  center,  at 
radial  distances  from  the  center  that  are  inside  the  granule,  at  the  granule  surface  and  at  radial 
distances  that  are  outside  the  granule. 

The  model  developed  by  Thompson  and  Gerstman  [13,14]  makes  several  simplifying 
assumptions  about  material  properties  and  about  the  nature  of  tissue-laser  interaction  in  the 
retina.  First,  it  assumes  melanosomes  are  homogenous  spheres  randomly  distributed  in  an 
infinite  volume  of  liquid  water.  Second,  it  assumes  absorption  of  laser  energy  takes  place  only 
within  melanin  granules  and  that  this  energy  absorption  is  homogeneous.  Third,  the  melanin 
granules  also  have  the  thermal  properties  of  water.  This  assumption  can  be  removed  in  favor  of 
a  more  accurate  solution  [9],  but  then  thennal  superposition  would  no  longer  be  valid. 

The  governing  equation  for  heat  transfer  from  a  single  granule  with  spherical  symmetry  is  shown 
in  Equation  (2.1). 


1  ^  on  i  Ao(r,t)  = 

r  dr 2  k  a  dt 


(2.1) 
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The  heat  equation  is  expressed  in  terms  of  r,  the  scalar  distance  from  the  granule  center  [cm],  T 
is  temperature  [°C],  t  is  time  [sec],  a  is  the  thermal  diffusivity  of  the  medium  [cm  /sec],  k  is  the 
thermal  conductivity  [J/cm  °C  sec],  and  A0  is  a  source  term  defined  in  Equations  (2.2)  and  (2.3). 


A0  =  — 

°  4ax 


1  - 


P  L 


2a„a2 


(l  -  exp(-2aam)  (1  +  2ama 


))] 


(2.2) 


A0 


if  0  >  r  >  a;  0  >  t  >  Tp 
otherwise 


(2.3) 


2 

The  variable  I0  is  the  average  retinal  fluence  [J/cm  ],  a  is  the  granular  radius  [pm],  tp  is  the  pulse 
length  [sec],  and  am  is  the  wavelength-dependent  melanosome  absorption  coefficient  [cm1]. 


2.2  FORTRAN  Code 

The  FORTRAN  implementation  of  the  Thompson-Gerstman  model  provided  computational 
analysis  of  damage  given  a  set  of  conditions  and  parameters  that  could  be  approximated  from 
data  gathered  experimentally.  A  rectangular  coordinate  system  was  used  to  define  the  locations 
of  melanin  granules  and  to  reference  the  spatial  points  at  which  temperature  rises  were 
computed.  The  origin  of  the  coordinate  system  is  defined  as  the  center  of  the  rectangular  volume 
shown  in  Figure  1. 
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Figure  1:  A  Profile  of  the  Distribution  of  Melanin  Granules  About  the  RPE  Layer 
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The  structure  of  the  FORTRAN  code  can  be  described  in  three  steps: 


1)  Input  Parameters:  The  initial  conditions  needed  for  the  temperature 
calculations,  as  well  as  a  set  of  example  values,  are  shown  in  Table  1. 
Values  related  to  simulated  domain  and  physical  parameters  are  used 
to  calculate  other  variables  needed  for  initialization.  These 
parameters  are  stored  in  “Config  Files”  shown  in  Figure  2. 

2)  Thermal  Solution:  A  temperature  value  is  calculated  for  each  grid 
point  the  user  specifies  over  an  interval  of  time  such  that  gradients  in 
thermal  diffusion  have  significantly  decreased.  The  time -temperature 
values  are  obtained  by  a  superposition  of  the  thermal  energy 
contributed  by  each  granule.  The  position,  time,  and  temperature  data 
is  then  saved  in  a  text  file  shown  in  Figure  2. 

3)  Damage:  Since  the  aim  of  the  model  is  to  determine  if  damage  has 
occurred  under  laser  radiation  exposure,  time -temperature  data  alone 
are  insufficient.  The  Arrhenius  damage  integral  was  selected  for  its 
usefulness  in  modeling  thermal  damage  to  biological  tissue.  The 
damage  data  are  saved  as  a  “Raw  Damage  File”  shown  in  Figure  2. 


Figure  2:  Flow  Diagram  of  the  FORTRAN  Thompson-Gerstman  Model 


Table  1:  Initial  Values  Needed  in  FORTRAN  Version  of  Thompson-Gerstman  Model 


Description 

Value 

X  domain 

[-20,20]  pm 

Y  domain 

[-20,20]  pm 

Z  domain 

[-7,7]  pm 

Number  of  melanin  granules 

186 

Simulation  time 

[0,2]  ms 

Comeal  fluence 

3  x  10'5  J/cm2 

Pulse  duration 

3  x  10'6  s 

Thermal  diffusivity  of  the  medium 

0.00133  cnfVs 

Thermal  conductivity  of  the  medium 

0.00556  J/cm-C-sec 

Beam  diameter  at  cornea 

0.7  cm 
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2.3  Source  Types 

The  type  of  source  exposure  can  also  be  designated  as  one  of  three  options:  gaussian,  top-hat,  or 
annular.  The  gaussian  profile  is  most  commonly  used  as  the  assumed  profile  of  most  laser 
exposures  and  represents  an  ideal  beam.  A  top-hat  profile  has  a  single  region  of  high  intensity 
surrounded  by  a  region  of  zero  intensity  without  the  smooth  transition  seen  in  the  gaussian 
profile.  An  annular  beam  resembles  a  ring  or  donut-shaped  beam  similar  to  the  top-hat  profile 
with  a  small  region  of  zero  intensity  in  the  center  of  the  beam  profile. 


Examples  of  thermal  results  using  the  parameters  in  Table  1  are  shown  in  Figure  3.  The  false 
color  images  show  temperature  rise  in  an  x-y  slice  at  the  beam  focus. 


,70  20 

60  15 

50„  10 
|  5 
40S 
u>  o 
30< 


-20-15-10  -5  0  5  10  15  20 
X-Axis  [pm] 


-20-15-10  -5  0  5  10  15  20 
X-Axis  [pm] 


20 

10 

0 


-5 

-10 

-15 

-20 


Figure  3:  Thermal  Profiles  Resulting  from  (a)  Gaussian  Beam,  (b)  Top-hat  Beam  and 

(c)  Annular  Beam 


2.4  Damage 

While  it  is  true  that  tissue  heated  to  the  vaporization  point  of  water  will  likely  cause  damage,  it  is 
also  possible  for  damage  to  occur  at  lower  temperatures.  Experiments  indicate  that  temperature 
elevations  of  just  Iff  C  can  cause  cellular  damage  if  maintained  for  a  long  enough  time.  The 
Arrhenius  damage  integral  approach  [15]  assumes  damage  occurs  by  an  active  rate  process 
described  by 

[l(r)=A  (2.4) 

The  variable  Q  represents  the  damage  level,  E  is  the  activation  energy  [energy  per  mole],  R  is 
the  gas  constant,  T  is  the  temperature,  and  tj  through  tf  represents  the  lower  and  upper  limit  of 
time  evaluated.  The  constants 

24  =  1.3x10"  [- 

.5. 

and 

E  =  1.5  x  105  [— 1 

LmoiJ 
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were  found  [16]  to  give  a  simple  value  for  fl  to  determine  if  damage  had  occurred. 

(damage  fl  >  1 

[no  damage  H  <  1 

Using  the  time-temperature  data  of  the  exposure  simulation,  the  Arrhenius  integral  is  solved 
numerically  by  Simpson’s  rule.  An  example  of  a  damage  spot  calculated  using  this  method  is 
shown  in  Figure  4.  The  damage  spot  was  calculated  by  the  Arrhenius  integral  for  exposure  to  a 
gaussian  beam  with  input  parameters  given  in  Table  1. 


X-Axis  [pm] 

Figure  4:  Damage  Spot  from  Gaussian  Beam 


3  IMPROVEMENT  OF  THOMPSON-GERSTMAN  MODEL 
3.1  C++  Implementation 

Several  improvements  were  suggested  to  the  original  FORTRAN  implementation  of  the 
Thompson-Gerstman  model.  First,  the  method  needed  to  be  converted  to  a  language  which 
allowed  for  added  usability  and  easy  integration  into  other  programs.  Second,  the  temperature 
data  should  be  passed  directly  to  the  damage  integral  instead  of  reading  and  writing  time- 
temperature  data  to  file  between  each  step.  Third,  the  user  should  be  able  to  provide  a  range  of 
values  as  input  parameters  instead  of  requiring  modifying  the  input  file  for  each  run.  Finally,  an 
option  of  multiple-pulse  exposures  could  be  added  to  expand  the  application  of  this  model 
beyond  the  single-pulse  realm.  With  these  criteria,  the  C++  language  was  selected  to  build  the 
new  implementation  of  the  model. 
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Figure  5:  Absolute  Temperature  Difference 


A  validation  of  the  C++  implementation  was  performed  by  direct  comparison  to  the  thermal 
results  given  by  the  FORTRAN  code.  This  thermal  results  comparison,  shown  in  Figure  5, 
produced  maximum  differences  in  temperature  around  0.08  °C.  This  gives  a  maximum  relative 
percent  error  of  0.04%.  This  difference  is  likely  due  to  minor  differences  between  how  the  two 
languages  handle  float  point  precision. 


3.2  C++  Code  Structure 


>  The  C++  version,  through  use  of  modular 
design,  was  able  to  significantly  streamline  the 
process  of  generating  damage  files.  Many  of 
these  changes  reflected  the  intended  usage  and 
computing  power  available  for  the  new 
implementation. 

>  This  version  of  the  Thompson-Gerstman  model 
was  intended  to  be  a  modular  tool  fit  for 
integration  into  other  computational  models. 
This  adds  usability  to  the  standalone  code 
which  can  now  be  ported  as  a  small  function 
library.  The  new  top-down  approach  of  the 
structure  is  shown  in  Figure  6. 


Figure  6:  Top  Down  Overview  of  C++  Version  of  Thompson-Gerstman  Model 
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3.3  Variable  Input  Parameters 

The  new  implementation  has  the  ability  to  iterate  the  thermal  solver  algorithm  over  any  number 
of  different  parameters  such  as  comeal  fluence,  laser  power,  number  of  points,  and  other  input 
parameters  of  significance.  This  modularized  code  might  be  used  to  optimize  input  parameters 
in  the  model  to  produce  effects  which  match  experimental  values.  Figure  7  shows  the  change  in 
temperature  in  the  medium  at  the  beam  focus  shortly  after  a  single  pulse.  A  range  of  corneal 
fluence  values  between  3.0  -  4.5  pJ/cnr  was  provided  while  holding  other  variables  constant. 
As  expected,  a  linear  relationship  between  comeal  fluence  and  temperature  exists  directly  after 
the  pulse. 
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Figure  7:  Temperature  Rise  at  Beam  Focus 


3.4  Multi-Pulse  Capability 

While  meaningful  analysis  can  be  perfonned  by  varying  parameters  for  a  single-pulse  exposure, 
multiple-pulse  trains  present  several  additional  parameters  to  compare  damage  thresholds.  Due 
to  the  unique  way  the  Thompson-Gerstman  model  handles  temperature  changes  from  a  single 
pulse,  implementing  a  multi-pulse  function  simply  requires  additional  iterations  of  the 
superposition  temperature  solver.  This  allows  for  variance  in  number  of  pulses  as  well  as  for  the 
time  between  pulses. 
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Figure  8:  Multi-Pulse  Validation 

Validation  of  this  method  was  performed  by  direct  comparison  to  another  multiple-pulse 
computational  solver,  which  has  been  independently  verified.  Figure  8  shows  the  direct 
comparison  of  the  multi-pulse  Thompson-Gerstman  algorithm  (green)  next  to  the  verified 
TempBuilder  (red)  thermal  profile.  These  results  show  good  agreement,  which  supports  the 
validity  of  the  new  C++  Thompson-Gerstman  code  for  multiple-pulse  exposures. 


4  DISCUSSION 

With  the  growing  number  of  laser  systems  being  used  by  both  military  and  civilians,  laser 
radiation  damage  to  the  retina  is  a  growing  concern.  Damage  to  tissue  and  loss  of  vision  can 
occur  under  a  variety  of  beam  conditions,  wavelengths,  and  exposure  times.  Developing  a  set  of 
laser  safety  standards  requires  testing  exposure  conditions  to  define  safe  limits  of  laser  use. 
Simulated  modeling  of  the  parameter  ranges  where  exposure  conditions  cause  photothermal 
damage  will  better  help  researchers  predict  damage. 


The  Thompson-Gerstman  model  is  able  to  give  an  estimate  of  whether  damage  has  occurred  by 
making  several  assumptions  that  are  not  physically  valid,  but  may  provide  model  results  that  are 
a  reasonable  match  to  actual  biological  damage  data.  First,  the  model  is  not  given  a  thermal 
boundary  condition,  so  the  exact  rate  of  thermal  energy  transmitted  to  surrounding  tissue  is 
neglected.  For  exposures  greater  than  1  second,  or  long  simulation  times  where  the  diffusion 
term  dominates  the  heat  equation,  this  assumption  would  likely  cause  the  model’s  temperature  to 
be  lower  than  expected.  Additionally,  since  photomechanical  and  photochemical  damage  are  not 
included  in  the  model,  short  pulses  and  temperatures  at  or  near  vaporization  should  be  ignored  as 
non-valid.  However,  since  photothermal  damage  covers  a  wide  range  of  possible  conditions, 
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other  damage  mechanisms  are  not  needed  as  long  as  the  model  is  used  in  a  valid  laser  pulse 
range. 

The  positive  aspects  of  this  model  lie  mostly  in  its  numerical  simplicity  and  speed  of  simulation. 
Since  computing  temperature  is  simplified  using  superposition  of  analytic  solutions,  long  lists  of 
simultaneous  equations  or  numerical  methods  of  solving  partial  differential  equations  are 
unnecessary.  As  a  result,  the  computational  effort  is  reduced  to  several  loops  of  simple 
arithmetic.  This  simplifies  what  could  be  an  exhaustive  parameter  study  into  a  series  of  non- 
intensive  computational  steps. 


5  CONCLUSIONS 

The  work  which  began  with  the  FORTRAN  implementation  of  the  Thompson-Gerstman  model 
has  been  expanded  to  a  C++  modular  model  that  possesses  a  function  library  capable  of 
interfacing  with  other  laser-tissue  damage  models.  This  modular  model  has  improved  the 
usability  of  the  method  for  conducting  damage  threshold  research.  The  granular  model  described 
here  and  detailed  in  their  paper  [14]  has  proven  to  be  a  fast  and  efficient  method  at  predicting 
retinal  damage  using  a  wide  range  of  input  parameters.  Additional  developments  of  multi -pulse 
functionality  and  solving  exposures  over  large  ranges  of  input  parameters  greatly  increase  the 
application  and  research  potential  of  this  method.  Direct  comparison  validation  has  been 
performed  on  the  results  of  the  FORTRAN  and  C++  codes,  here  showing  good  agreement 
between  the  two  computational  codes.  Results  from  this  model  could  be  compared  to  known 
biological  results,  which  would  possess  insight  into  thresholds  in  regimes  where  no  laser  damage 
studies  have  been  conducted.  This  tool  would  have  the  added  benefit  of  reducing  the  complexity 
and  sensitivity  of  attaining  results  when  depending  on  animal  studies  for  determining  damage 
thresholds. 
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7  APPENDIX:  THOMPSON  GERSTMAN  MODEL  C++ 


7.1  Top  Level  Function 

#include<math . h> 

#include<ctime> 

#include  <boost/test/unit_test . hpp> 
#include<iostream> 

#include<f stream> 

#include<vector> 

#include<string> 
using  namespace  std; 

#include  "TGlib.cpp" 

#include  <boost/property_tree/ptree . hpp> 

# include  <boost/property_tree/ ini_parser . hpp> 
#include<omp . h> 

int  main(int  argc,  char*  arg [ ] ) { 

BigTherm  Therm (arg [ 1 ]) ; 

Therm. TakataDamage ( ) ; 
return  0; 

} 
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7.2  Header  File 


#include<vector> 

#include  <boost/property_tree/ptree . hpp> 

#include<string> 

using  namespace  std; 

class  granule  { 
public : 


// -  These  parameters  are  given  in  . inpl  files - 

double  a;  //Granular  radius 

// -  End  parameters  given  by  .inpl  files - 

// -  These  parameters  are  givin  in  . inp2  files  - 

double  pulsedur;  //pulse  duration  in  seconds 

double  alpha;  //  Thermal  Diffusivity  of  medium  in  cm**2/sec 

double  cond;  //  Thermal  conductivity  of  medium  in  j/cm-Csec 

// -  End  Parameters  from  . inp2  files  - 


//used  in  (orig, in, surf , out, vanilla) temp  functions,  calculated  in 
//bigthrem  from  other  input  parameters 

double  AO; 

double  OutTemp (double  r,  double  t) ; 
double  Surf Temp (double  t)  ; 
double  InTemp (double  r,  double  t)  ; 
double  OrigTemp (double  t)  ; 
double  Temp (double  r,  double  t) ; 
double  MelTemp (double  r,  double  t)  ; 
double  tmin,  PulseSep; 
int  PulseNum; 
vector<double>  PulseVec; 

}  ; 


class  BigTherm  { 
public : 

BigTherm (char*  conf igf ilename) ; 
vector<int>  elapsed,  tarray,  guess; 
int  k,  s,  improf,  n,TempFlag,  SeedFlag; 
int  xnum,  ynum,  znum,  tnum,  melnum,  rnum; 
vector<int>  test; 

double  BoxLength,  BoxDepth,  MelDens; 
double  xmin,  xmax,  xsize; 
double  ymin,  ymax,  ysize; 
double  zmin,  zmax,  zsize; 
double  tmin,  tmax,  tsize; 
vector<  vector<  double  >  >  Pos; 
vector<  double  >  weight,  PulseVec; 
vector<  vector<  double  >  >  VecTemp; 
vector<  double  >  rsize; 
double  dx; 

double  rani,  ran2,  ran3; 
double  irr,Tbody,  absorb; 
double  Temp,  x,  y,  z,  t,  r; 
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double  sigma,  corspot.  Focus; 
double  spotsize, MelTemp, trans, dobs; 
double  Qtot,  Qlost,  Qkept; 
granule  granny; 

boost: : property_tree : :ptree  configfile; 
string  MelPlacementFile; 
string  TempOutputFile; 
string  DamageOutputFile; 

double  BigThermTempAt (double  x,  double  y,  double  z,  double  t) ; 

//  this  function  is  to  build  the  function  again  if  you  change  any 
//  variables  that  affect  weight  or  vectemp 
void  Reconstruct (void) ; 
void  BigTherm2 (void) ; 

void  BigTherm2 (double  outputxmin,  double  outputxmax, int  outxnum, 
double  outputymin,  double  outputymax,  int  outynum,  double  outputzmin,  double 
outputzmax,  int  outznum) ; 

int  PulseNum; 

double  PulseSep; 

void  TakataDamage (void) ; 

int  TakataDamageAt (double  x,  double  y,  double  z); 
int  Damagelnfo; 
int  numthreads; 

}  ; 

void  takatadamage (char*  conf igf ilename) ; 
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7.3  Library  File 


#include<math . h> 

#include"TGlib . h" 

#def ine  JJSE_MATH_DEFINES 
#include<limits> 

#include<vector> 

#include<f stream> 

#include  <boost/property_tree/ptree . hpp> 
#include  <boost/property  tree/ini_parser . hpp> 
#include<string> 

#include<iostream> 

#include" omp . h" 

#include  <time.h> 
using  namespace  std; 

double  Gmma (double  aa,  double  bb,  double  cc) { 


if (bb  ==  0.0)  { 

// - gigamma  function - 

if (cc  >  30.0) 

cc  =  30.0; 

// - gammaser  function - 

// - modified (un-normalized)  version  of  GSER  from  numerical  recipes  - 

/* 


Modification  of  GSER  is  ignoring  (removing)  the  "gin"  variable,  and  the 
function  that  calculates  it  gammln(a) .  this  is  used  in  the  sum  evaluating  the 
exponential  for  del  <  sum  *EPS,  as  subtraction  by  the  natural  log  of  the 
complete  gamma  function  inside  the  exponential  acts  the  same  as  dividing  by 
the  same  value. 

*/ 

const  int  ITMAX  =  100; 

const  double  EPS  =  numeric  limits<double> :: epsilon () ; 
double  sum,  del,ap; 
if (cc  <=  0 . 0 ) { 

if (cc<0 . 0 ) { 

return  -1; 

//fail  in  a  mysterious  way 

} 

else 

return  0.0; 

} 

else  { 

ap  =  aa; 

del  =  sum  =  1.0/aa; 
for (int  n  =  0;  n  CITMAX;  n++) { 
ap  +=  1.0; 
del  *=  cc/ap; 
sum  +=  del; 

if(fabs(del)  <  f abs ( sum) *EPS ) { 

return  sum*exp (-cc+aa*log (cc)  )  ; 

} 

} 

} 

} 

} 

double  LagL (double  a,  double  b,  double  c) { 
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double  d,  gml,  gm2.  Gamma,  answer; 
d  =  -c; 
gml  =  -0.5; 
gm2  =  0.0; 

Gamma  =  Gmma (gml , gm2 , d) ; 

answer  =  -  (Gamma  *  sqrt (d)  *  (1.0  -  2*c/3)  +  (2*exp (c) ) /3) /M_PI; 

return  answer; 

} 

//calculates  the  temperature  outside  the  granule 
double  granule :: OutTemp (double  r,  double  t) { 
double  Outterm[4],  gm[2],  lg[2]; 
double  OutTheta; 
gm [ 0 ]  =  0.5; 
gm [ 1 ]  =  0.0; 
lg [ 0 ]  =  1.5; 
lg [ 1 ]  =  -0.5; 

double  gamma  =  Gmma (gm[0] , gm [ 1 ] , ( ( (-a  +  r) * (-a  +  r) ) / (4 . 0*alpha*t) ) ) ; 

Outterm[0]  =  ( (a*t* ( (sqrt (M_PI) * ( (-a  +  r)*(-a  +  r)  + 

2 . 0*alpha*t) ) / (2 . 0*alpha*t)  -  (2.0*(-a  +  r ) * ( 0 . 5*exp ( ( - (  ( -a  +  r)*(-a 

+r) ) / (4 . 0*alpha*t) ) )  +  (sqrt (((-a  +  r)*(-a  +r) ) / (alpha*t) ) * (  (  (-a  +  r)*(-a 

+r)  )  +  2 . 0 * alpha* t) *  gamma) / (4 . 0*  (  (-a  +  r)  *  (- 

a+r)  )  )  )  )  / (sqrt (alpha) *sqrt (t) ) ) ) / (sqrt (M_PI) ) ) ; 

Outterm[l]  =  (a*a*a  -  3.0*a*a*r  +  3.0*a*r*r  -  (r*r*r)  +  6 . 0*a*alpha*t  - 

6 . 0*alpha*r*t  + 

6.0* (alpha* sqrt (alpha) ) *sqrt (M_PI ) * ( t*sqrt ( t ) ) *LagL (lg[0],lg[l],-((-a  +  r) * (- 

a  +  r) ) / (4 . 0*alpha*t) ) ) / (6 . 0*alpha) ; 

Outterm[2]  =  (-a*a*a  -  3.0*a*a*r  -  3.0*a*r*r  -  r*r*r  -  6 . 0*a*alpha*t  - 

6 . 0*alpha*r*t  + 

6.0* (alpha* sqrt (alpha) ) *sqrt (M_PI ) * ( t*sqrt ( t ) ) *LagL ( lg [ 0 ] , lg [ 1 ]  ,  (-(a  +  r ) * (a 

+r) ) / (4 . 0*alpha*t) ) ) / (6 . 0* alpha) ; 

Outterm[3]  =  (a*t* ( (sqrt (M_PI)  *  ((a  +  r) * (a  +  r)  + 
2 . 0*alpha*t) ) / (2 . 0*alpha*t)  -  (2.0*(a  +  r) *  (0 . 5*exp  (- ( (a  +  r) *  (a  + 
r) ) / (4 . 0*alpha*t) )  +  (sqrt (((a  +  r) * (a  +  r)  ) / (alpha*t) ) * ( (a  +  r) * (a  +  r)  + 
2 . 0*alpha*t) *Gmma (gm[0] ,  gm[l],  ((a  +  r) * (a  +  r) ) / (4 . 0*alpha*t) ) ) / (4 . 0* (a  + 
r)  *  (a  +  r)  )  )  )  /  (sqrt  (alpha)  *sqrt  (t)  )  )  )  /  (sqrt  (M_PI)  )  ; 

OutTheta  =  alpha*A0* (Outterm[0]  -  Outterm[l]  +  Outterm[2]  + 
Outterm[3] ) / (2 . 0*cond)  ; 
return  OutTheta/r; 

} 

//  calculates  temperature  on  the  surface  of  the  granual 

double  granule :: Surf Temp (double  t) { 

double  Surffnc[4]; 

double  SurfTheta,  gm[2],  lg[2]; 

gm [ 0 ]  =  0.5; 

gm [ 1 ]  =  0.0; 

lg [ 0 ]  =  1.5; 

lg [ 1 ]  =  -0.5; 

Surffnc[0]  =  t; 

Surffnc[l]  =  t*(sqrt(M  PI)*(2.0*a*a  +  alpha*t) / (alpha*t)  -  4 . 0*a*  (0 . 5*exp (- 

(a*a) / (alpha*t) )  +  sqrt ( (a*a) / (alpha*t) ) * (2 . 0*a*a  + 

alpha* t) *Gmma (gm[0] , gm [ 1 ] , (a* a) / ( alpha* t) ) / (4.0*a*a) ) / (sqrt (alpha) *sqrt (t) ) ) / 
(sqrt (M_PI ) ) ; 

Surffnc[2]  =  4 . 0* (t*sqrt (t) ) / (3 . 0*sqrt (M_PI) ) ; 
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Surffnc[3]  =  -4 . 0* (a*a*a) / (3 . 0*alpha*sqrt (alpha) )  -  2 . 0*a*t/sqrt (alpha)  + 

sqrt (M_PI)  *t*sqrt (t) *LagL (lg [0] , lg [1] , - ( (a*a) / ( alpha* t) ) ) ; 

Surffnc[4]  =  a*alpha*A0* (Surffnc  [0]  +  Surf fnc [ 1 ] ) / (2 . 0*cond)  + 

alpha*sqrt (alpha) *A0* (Surf fnc [3]  -  Surf fnc [2 ] ) / (2 . 0*cond) ; 
return  Surf fnc [ 4 ] /a; 

} 

//  Calculates  temperature  inside  the  granual 

double  granule :: InTemp (double  r,  double  t) { 

double  Interm [5],  gm[2],  lg[2]; 

double  Intheta; 

gm [ 0 ]  =  0.5; 

gm [ 1 ]  =  0.0; 

lg [ 0 ]  =  1.5; 

lg [ 1 ]  =  -0.5; 

lnterm[0]  =  r*t; 

Interm[l]  =  (-a*a*a  +  3.0*a*a*r  -  3.0*a*r*r  +  r*r*r  -  6 . 0 *a*alpha*t  + 

6 . 0*alpha*r*t  +  6 . 0*alpha*sqrt (alpha) *sqrt (M  PI) *t*sqrt (t) *LagL (lg [0] , lg [1] , - 

((a  -  r) * (a  -  r) ) / (4 . 0*alpha*t) ) ) / (12 . 0*alpha) ; 

Interm[2]  =  (-a*a*a  -  3.0*a*a*r  -  3.0*a*r*r  -  r*r*r  -  6 . 0*a*alpha*t 

6 . 0*alpha*r*t  +  6 . 0*alpha*sqrt (alpha) *sqrt (M  PI) *t*sqrt (t) *LagL (lg [0] , lg  [1] , - 

((a  +  r) * (a  +  r)  ) /  (4 . 0*alpha*t) ) ) / (12 . 0*alpha)  ; 

Interm[3]  =  a*t*(sqrt(M  PI)*(((a  -  r) * (a  -  r) )  +  2 . 0*alpha*t) / (2 . 0*alpha*t)  - 

2.0*(a  -  r) * (0 . 5*exp (- ( (a  -  r) * (a  -  r) ) / (4 . 0*alpha*t) )  +  (a  - 

r) /sqrt ( alpha* t) *((a  -  r)*(a  -  r)  +  2. 0*alpha*t) *Gmma (gm[0] , gm [ 1 ]  ,  ( (a  -  r) * (a 
r)  ) / (4 . 0*alpha*t) ) / (4 . 0* (  (a  -  r) *  (a 

r)  )  )  ) /sqrt ( alpha* t) ) / (2 . 0*sqrt (M_PI ) ) ; 

Interm[4]  =  a*t*(sqrt(M  PI)*(((a  +  r) * (a  +  r)  )  +  2 . 0*alpha*t) / (2 . 0*alpha*t)  - 

2.0*(a  +  r) * (0 . 5*exp (- ( (a  +  r) * (a  +  r) ) / (4 . 0*alpha*t) )  +  (a  + 

r)  /sqrt ( alpha* t) *(((a  +  r)*(a  +  r)  )  +  2. 0 * alpha* t) *Gmma (gm[0] , gm  [  1 ] ,  (  (a  + 

r) *  (a  +  r)  ) / (4 . 0*alpha*t) ) / (4 . 0*  (  (a  +  r) * (a  + 

r)  )  )  )  /sqrt ( alpha* t) ) / (2 . 0*sqrt  (M_PI ) ) ; 

Intheta  =  alpha*A0* (lnterm[0]  -  Interm[l]  +  Interm[2]  -  Interm[3]  + 

Interm[4] ) /cond; 

return  Intheta/r; 

} 

double  granule :: OrigTemp (double  t) { 
double  Origfnc[3],  gm[2],  lg[2]; 
gm [ 0 ]  =  0.5; 
gm [ 1 ]  =  0.0; 
lg [ 0 ]  =  1.5; 
lg [ 1 ]  =  -0.5; 

Origfnc[0]  =  t; 

Origfncfl]  =  t*(sqrt(M  PI)*(a*a  +  2 . 0*alpha*t) / (2 . 0*alpha*t) 

2 . 0*a* (0 . 5*exp (-a*a/ (4 . 0*alpha*t) )  +  sqrt (a*a/ (alpha*t) ) * (a*a  + 

2 . 0 * alpha* t) *Gmma (gm[0] , gm [ 1 ] , a* a/ (4 . 0 * alpha* t) ) / (4.0*a*a) ) / (sqrt (alpha) *sqrt 
(t) ) ) /sqrt (M_PI ) ; 

Origfnc[2]  =  sqrt (t) *(- (a*sqrt (M_PI) / (sqrt (alpha) *sqrt (t) ) ) 

sqrt (a* a/ ( alpha* t) ) *Gmma (-gm[0] , gm [ 1 ] , a* a/ (4 . 0 * alpha* t) ) /2 . 0) /sqrt (M_PI) ; 
return  alpha*A0* (Origfnc [0]  -  Origfnc[l]  -  (a*Origfnc [2 ]) /sqrt (alpha) ) /cond; 

} 

double  granule :: Temp (double  r,  double  t) { 
if (t  ==  0.0) 

return  0.0; 
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else  { 


if (r  ==  0.0) 

return  OrigTemp(t); 

else  { 

if ( (r>0 . 0 ) && (r<a) ) 

return  InTemp(r,t); 

else  { 

if  (r  ==  a) 

return  SurfTemp(t); 

else  { 

if (r>a) 

return  OutTemp (r, t) ; 

} 

} 

} 

} 

} 

double  granule :: MelTemp (double  r,  double  t) { 
double  times; 
double  ReturnTemp; 

for(int  i  =0;  i<PulseVec . size ( ) ;  i++) { 
times  =  t  -  PulseVecfi]; 

/* 

This  part  here  is  basically  the  original  MelTemp  Function  from  FORTRAN,  here 
it  is  looped  over  the  number  of  pulses  that  are  given 
*/ 

if (  times  <  0 . 0 ) { 

ReturnTemp  +=  0.0; 

} 

else  { 

if (times  >  pulsedur) 

ReturnTemp  +=  Temp (r, times)  -  Temp (r, times-pulsedur)  ; 

else 

ReturnTemp  +=  Temp ( r, times ) ; 

} 

} 

return  ReturnTemp; 

} 

BigTherm: : BigTherm (char*  conf igf ilename) { 
boost: : property_tree : :ptree  conf igf ile; 

boost : : property_tree : : ini_parser : : read  ini (conf igf ilename,  configfile) ; 
BoxLength=conf igf ile . get<double> ( "Bigtherm. BoxLength" ) ; 

BoxDepth=conf igf ile . get<double> ( "Bigtherm. BoxDepth" )  ; 

granny . a=conf igf ile . get<double> ( "Bigtherm. GranularRadius" ) ; 

xnum=conf igf ile . get<int> ( "Bigtherm. xnum" ) ; 

xmin=conf igf ile . get<double> ( "Bigtherm. xmin" ) ; 

xmax=conf igf ile . get<double> ( "Bigtherm. xmax" ) ; 

ynum=conf igf ile . get<int> ( "Bigtherm. ynum" ) ; 

ymin=conf igf ile . get<double> ( "Bigtherm. ymin" ) ; 

ymax=conf igf ile . get<double> ( "Bigtherm. ymax" ) ; 

znum=conf igf ile . get<int> ( "Bigtherm. znum" ) ; 

zmin=conf igf ile . get<double> ( "Bigtherm. zmin" ) ; 

zmax=conf igf ile . get<double> ( "Bigtherm. zmax" ) ; 

tnum=conf igf ile . get<int> ( "Bigtherm. tnum" ) ; 

tmin=conf igf ile . get<double> ( "Bigtherm. tmin" ) ; 

19 

Distribution  A:  Approved  for  public  release;  distribution  unlimited.  PA  Case  No:  TSRL-PA-20 15-0013. 


tmax=conf igf ile . get<double> ( "Bigtherm. tmax" ) ; 

MelDens=conf igf ile . get<double> ( "Bigtherm. Mel Dens" ) ; 
rnum=conf igf ile . get<int> ( "Bigtherm. rnum" ) ; 
absorb=conf igf ile . get<double> ( "Bigtherm. absorb" ) ; 
irr=conf igf ile . get<double> ( "Bigtherm. irr" ) ; 

granny . pulsedur=conf igf ile . get<double> ( "Bigtherm. pul sedur" ) ; 
granny . alpha=conf igf ile . get<double> ( "Bigtherm. alpha" ) ; 
granny . cond=conf igf ile . get<double> ( "Bigtherm. cond" ) ; 
corspot=conf igf ile . get<double> ( "Bigtherm. cor spot" ) ; 
spotsize=conf igf ile . get<double> ( "Bigtherm. spot size" ) ; 
trans=conf igf ile . get<double> ( "Bigtherm. trans" )  ; 
improf =conf igf ile . get<int> ( "Bigtherm. improf " ) ; 

SeedFlag=conf igf ile . get<int> ( "Bigtherm. SeedFlag" ) ; 

PulseNum=conf igf ile . get<int> ( "Bigtherm. PulseNum" ) ; 

PulseSep=conf igf ile . get<double> ( "Bigtherm. PulseSep" ) ; 

PulseVec . push  back(tmin); 
for(int  i  =1;  i<PulseNum;  i++) { 

PulseVec . push  back ( PulseVec [ i-1 ]  +  granny . pulsedur  +  PulseSep); 

} 

granny . PulseVec  =  PulseVec; 

// - 

TempOutputFile=conf igf ile . get<string> ( "Bigtherm. TempOutputFile" ) ; 
DamageOutputFile=conf igf ile . get<string> ( "Bigtherm. DamageOutputFile" ) ; 
Damage Inf o=conf igf ile . get<int> ( "Bigtherm. Damage Info" ) ; 
numthreads=conf igf ile . get<int> ( "BigTherm. numthreads" , 1 ) ; 
if (numthreads<l ) 

numthreads  =  1; 

////////////////////////////////////////////////////////////////// 

/ / end  loading 

////////////////////////////////////////////////// 

xsize  =  (xmax  -  xmin) /xnum; 

ysize  =  (ymax  -  ymin)/ynum; 

zsize  =  (zmax  -  zmin) /znum; 

tsize  =  (tmax  -  tmin)/tnum; 

melnum  =  0; 

melnum  =  conf igf ile . get<int> ( "Bigtherm. Melnum" ,  0); 
if (melnum  ==  0){ 

MelDens=conf igf ile . get<double> ( "Bigtherm. Mel Dens" )  ; 
melnum  =  MelDens  *  BoxLength*BoxLength  *  BoxDepth; 

} 

else 

MelDens  =  0; 

/* 

I  ■k-k-k-k'k-k-k'k'k-k'k-k'k-k-k-k'k-k'k-k'k-k-k-k'k-k'k-k'k-k-k'k'k-k'k-k'k-k-k-k'k-k-k'k'k-k'k-k-k'k'k'k-k'k-k'k'k'k-k'k-k'k-k-k-k'k-k-k-k'k-k 

I 

!  Determine  the  image  profile  on  the  retina. 

!  Use  image  profile  technique  determined  by 

!  the  value  of  improf: 


improf=0 

<-> 

Gaussian 

Image 

improf =3 

<-> 

Top  Hat 

Source 

improf=4 

<-> 

Annular 

Beam 

The  image  diameter  is  then  determined  by  the  value  of 
spotsize.  The  amount  of  focusing 

that  has  taken  place  is  a  function  of  the  retinal  image  diameter 
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!  as  well.  For  annular  beams,  which  assume  a  37  per  cent  central 

!  obscuration,  the  diameter  of  the  central  obscuration,  dobs,  is 

!  spotsize*DSQRT (0 . 37) . 

I 

I  k kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
*/ 


dobs  =  spotsize  *  sqrt(0.37); 
if ( improf  ==  0  ) 

Focus  =  2  *  (corspot*corspot)  /  (spotsize*spotsize) ; 
if  (improf  ==  3) 

Focus  =  (corspot*corspot)  /  (spotsize*spotsize) ; 
if ( improf  ==  4 ) 

Focus  =  (corspot*corspot)  /  (spotsize*spotsize) ; 
granny. AO  =  3  *  irr  *  trans  *  Focus  *  (1  -  (1  -  exp  (-2  * 

granny. a  *  absorb)*  (1  +  2  *  granny. a  *  absorb) )/ (2  *  granny . a*granny . a  * 
absorb*absorb)  )/  (4  *  granny. a  *  granny . pulsedur ) ; 

sigma  =  8.0  /  (spotsize*spotsize) ; 

if  (tmax  <  granny . pulsedur ) 

Qtot  =  granny. AO  *  4.0  *  M  PI  *  granny . a*granny . a*granny .  a  *  tmax  / 

3.0; 

else 

Qtot  =  granny. AO  *  4.0  *  M  PI  *  granny .  a*granny .  a*granny .  a  * 

granny . pulsedur  /  3.0; 

vector<double>  dummy; 

/ /  computing  weighting  of  the  f luence  for  the  placed  granule 
if(SeedFlag  ==  0){ 

MelPlacementFile=conf igf ile . get<std : : string> ( "Bigtherm. Mel Placement File 

")  ; 

f stream  posf ile (MelPlacementFile . c_str ()); 
if (pos file . good ( ) ) { 
double  in; 
string  throwaway; 
while ( Iposfile.eof  () )  { 
posf ile>>in; 
dummy. push  back(in); 
posf ile>>in; 
dummy. push  back(in); 

posf ile>>in; 
dummy. push  back(in); 

Pos. push  back(dummy); 

dummy . clear ( ) ; 

getline (posf ile, throwaway) ; 

} 

} 

else  { 

cout<<"Placement  File  not  found,  using  random  seed"<<endl; 

SeedFlag  =  1; 

} 

} 

double  thermcon; 
if (SeedFlag  ==  1){ 

srand (time (NULL) ) ; 

for (  int  j  =0;  j<melnum;  j++)  { 
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if  (j 


==  0) 

Pos . clear ( ) ; 
rani  =  (rand()  %  101)/100.0; 
ran2  =  (rand()  %  101J/100.0; 
ran3  =  (rand()  %  101J/100.0; 

thermcon  =  (rani  *  BoxLength)  -  (BoxLength/2 . 0 ) ; 
dummy. push  back (thermcon) ; 

thermcon  =  (ran2  *  BoxLength)  -  (BoxLength/2 . 0 ) ; 
dummy. push  back (thermcon)  ; 

thermcon  =  (ran3  *  BoxDepth)  -  (BoxDepth/2 . 0 ) ; 
dummy. push  back (thermcon) ; 

Pos. push  back(dummy); 
dummy . clear ( )  ; 

} 

} 

for (  int  j  =  0;  j<melnum;  j+L)  { 

r  =  sqrt (Pos  [  j ]  [0] *Pos  [  j ]  [0]  +  Pos  [  j ]  [1] *Pos  [  j ]  [1]  )  ; 
if  (improf  ==  0) 

weight. push  back (exp (-sigma  *  r*r) ) ; 
else  if  (improf  ==  3)  { 

if  (r  <=  spotsize/2 . 0 ) { 

weight. push  back (1.0); 

} 

else  { 

weight .  push_back (0.0) ; 

} 

} 

else  if  (improf  ==  4) { 

if (  r  >  spotsize/2 . 0 ) 

weight . push_back (0.0) ; 

else  { 

if (r  <  dobs/2 . 0 ) 

weight . push_back (0.0) ; 

else 

weight. push  back (1.0); 

} 

} 

} 

for (int  i  =  0;  i<=tnum;  i++) { 

t  =  tmin  +  i  *  tsize; 
rsize . push_back (0.0) ; 

if  (t  <=  0.00001) 

rsize[i]  =  (0 . 0005/rnum) ; 
else 

rsize[i]  =  (0.0005* (7.0  +  loglO (t) ) / rnum) ; 

for (int  j  =  0;  j<=rnum;  j++){ 
r  =  j  *  rsize  [i] ; 

dummy .push  back (granny . MelTemp (r, t) ) ; 
if(dummy[j]  +  Tbody  >  100) 

TempFlag  =  1; 


} 
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VecTemp.push  back(dummy); 
dummy . clear ( ) ; 

} 

//  -  Check  conservation  of  energy  - 

if  (tmax  <=  10.0E-6) 

dx  =  0.0005/rnum; 

else 

dx  =  0.0005*  (7.0  +  loglO (tmax) ) / rnum; 

Qkept  =  0.0; 

for(int  j  =0;  j<rnum;  j++){ 

Qkept  =  Qkept  +  VecTemp [tnum]  [j]  *  4.1868  *  4.0  *  M  PI  *  dx*dx*dx  * 

(3.0* (j+1) * (j+1)  -  3.0* ( j+1)  +  1.0)  /  3.0; 

} 

Qlost  =  Qtot  -  Qkept; 

if  (TempFlag  ==  1)  { 

cout<<endl ; 

cout<<"/ ///////////////////////////////////////////// /"«endl; 


cout«" - "<<endl  ; 

cout«"  Warning  .  .  .  100  degrees  exceeded.  "«endl; 

cout«" - "<<endl  ; 


cout<<"/ ///////////////////////////////////////////// /"«endl; 

} 


double  BigTherm: : BigThermTempAt (double  x,  double  y,  double  z,  double  t) { 
//WARNING  this  function  will  NOT  give  outputs  for  arbitrary  t,  simply  the 
//closest  t  to  the  tsize  given  away  from  the  min  t. 
if ( (t<tmin) | | (t>tmax) ) 
return  -1; 


int  1  =  (t-tmin) /tsize; 

double  dist,  ilook.  Temper,  SumTemp; 

SumTemp=0 . 0 ; 

//! -  Summing  effects  of  all  granules - 

omp_set  num  threads (numthreads ) ; 

#pragma  omp  parallel  for  reduction (+: SumTemp)  private (dist,  ilook.  Temper) 
for (int  m  =  0;  m<melnum;  m++) { 

dist  =  sqrt ( (Pos  [m]  [0] -x) * (Pos  [m]  [0] -x)  +  (Pos [m]  [1] -y) * (Pos  [m]  [1] -y)  + 

(Pos [m]  [2]-z)*(Pos[m]  [2] — z) ) ; 

/ / ! -  Individual  Temps  Interpolated  From  VecTemp  Table - 

ilook  =  dist/rsize  [1] ; 


if  (ilook  >  rnum-1)  { 

Temper=0 . 0 ; 

} 

else  { 

Temper  =  VecTemp [1] [ilook]  +  (dist  -  ilook  *  rsize[l]) 
(VecTemp[l] [ilook+1]  -  VecTemp[l] [ilook])  /  rsize[l]; 

} 

SumTemp  +=  weight [m]  *  Temper; 


} 

/* 
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Slight  issues  with  noise,  method  only  accurate  to  about  a  hundredth  of  a 
degree  anyway  so  throwing  out  temperature  rises  on  this  magnitude  cleans  up 
the 

profiles  a  bit. 

*/ 

if(SumTemp  <  le-6) 
return  0.0; 

return  SumTemp; 

} 

//This  function  is  a  copy  of  the  constructor  if  you  want  to  run  it  again  with 
the  //different  beam  parameters 

void  BigTherm: :ReConstruct (void) { 
double  thermcon; 
vector<double>  dummy; 

Pos . clear ( ) ; 

srand (time (NULL) )  ; 

for (  int  j  =  0;  j<melnum;  j++)  { 

if ( j==  0) 

Pos . clear ( ) ; 

rani  =  (rand()  %  101J/100.0; 

ran2  =  (rand()  %  101J/100.0; 

ran3  =  (rand()  %  101J/100.0; 

thermcon  =  (rani  *  BoxLength)  -  (BoxLength/2 . 0 ) ; 
dummy. push  back (thermcon)  ; 

thermcon  =  (ran2  *  BoxLength)  -  (BoxLength/2 . 0 ) ; 
dummy. push  back (thermcon) ; 

thermcon  =  (ran3  *  BoxDepth)  -  (BoxDepth/2 . 0 ) ; 
dummy. push  back (thermcon) ; 

Pos. push  back(dummy); 
dummy . clear ( )  ; 

} 

PulseVec . clear ( )  ; 

PulseVec . push  back(tmin); 
for(int  i  =1;  i<PulseNum;  i++) { 

PulseVec . push  back (PulseVec [ i —  1 ]  +  granny . pulsedur  +  PulseSep) ; 

} 

granny . PulseVec  =  PulseVec; 
weight . clear ( ) ; 

for (  int  j  =  0;  j<melnum;  j++)  { 

r  =  sqrt (Pos [ j ] [0] *Pos [ j ] [0]  +  Pos [ j ] [1] *Pos [ j ] [1]  )  ; 
if  (improf  ==  0) 

weight .push_back (exp (-sigma  *  r*r)  )  ; 
else  if  (improf  ==  3)  { 

if  (r  <=  spotsize/2 . 0 ) { 

weight. push  back (1.0); 

} 

else  { 

weight . push_back (0.0) ; 

} 

} 

else  if  (improf  ==  4) { 

if (  r  >  spotsize/2 . 0 ) 

weight . push_back (0.0) ; 
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else  { 

if (r  <  dobs/2 . 0 ) 

weight . push_back (0.0) ; 

else 

weight. push  back (1.0); 

} 

} 

} 

VecTemp . clear ( ) ; 
for(int  i  =  0;  i<=tnum;  i++) { 
t  =  tmin  +  i  *  tsize; 
rsize . push_back (0.0) ; 
if  (t  <=  0.00001) 

rsize[i]  =  (0 . 0005/rnum) ; 
else 

rsize[i]  =  (0.0005*  (7.0  +  loglO (t) ) / rnum) ; 
for(int  j  =  0;  j<=rnum;  j++){ 
r  =  j  *  rsize  [i] ; 

dummy .push  back (granny . MelTemp (r, t) ) ; 
if(dummy[j]  +  Tbody  >  100) 

TempFlag  =  1; 

} 

VecTemp. push  back(dummy); 
dummy . clear ( ) ; 


} 

//  -  Check  Conservation  of  Energy  - 

if  (tmax  <=  10.0E-6) 

dx  =  0. 0005/rnum; 

else 

dx  =  0.0005*  (7.0  +  loglO (tmax) ) /rnum; 

Qkept  =  0.0; 

for(int  j  =0;  j<rnum;  j++){ 

Qkept  =  Qkept  +  VecTemp [tnum]  [j]  *  4.1868  *  4.0  *  M  PI  *  dx*dx*dx  * 

(3.0* (j+1) *  (j+1)  -  3.0*  ( j+1)  +  1.0)  /  3.0; 

} 

Qlost  =  Qtot  -  Qkept; 

if  (TempFlag  ==  1)  { 

cout<<endl ; 

cout«"  //////////////////////////////////////////////////  "«endl  ; 


cout«" - "<<endl  ; 

cout«"  Warning  .  .  .  100  degrees  exceeded . "<<endl ; 

cout«" - "<<endl  ; 


cout«"  /  /////////////////////////////////////////////////  "«endl  ; 

} 

} 

void  BigTherm: :BigTherm2 (void) { 

/* 

This  will  run  the  equivalent  of  the  Original  bigtherm2,  and  write  it  to  the 
file  listed  for  output  in  the  config  file.  It  writes  the  entirety  of  the  box 
and  the  entirety  of  the  given  time. 

*/ 
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of stream  writef ile (TempOutputFile . c_str ( ) , ios : : out)  ; 
writef ile<<"t\t"<<"x\t"<<"y\t"<<"z\t"<<"Temp"<<"\n" ; 
for(int  i  =  0;  i<xnum+l;  i++) { 
x  =  xmin  +  i  *  xsize; 
for(int  ii  =  0;  ii<ynum+l;  ii++) { 
y  =  ymin  +  ii*ysize; 

for(int  iii  =  0;  iii<znum+l;  iii++) { 
z  =  zmin  +  iii*zsize; 
for(int  iv  =  0;  iv  <  tnum+1;  iv  ++) { 
t  =  tmin  +  iv  *  tsize; 

double  thing  =  BigThermTempAt (x, y, z, t) ; 
writef ile<<t<<" \t"<<x<<" \t"<<y<<" \t"<<z<<" \t"<<thing<<endl ; 


/* 

This  runs  BigTherm2  for  a  given  subset  of  the  box.  This  can  be  used  so  that 
one  does  not  output  or  calculate  an  unneeded  data,  i.e.  points  far  out  enough 
from  the  beam  that  there  is  no  significant  temperature  rise 
*/ 

void  BigTherm: : BigTherm2 (double  outxmin,  double  outxmax,  int  outxnum,  double 
outymin,  double  outymax,  int  outynum,  double  outzmin,  double  outzmax,  int 
outznum) { 

of stream  writef ile (TempOutputFile . c_str ( ) , ios : : out) ; 

writef  ile«"t\t"«" z\t"<<"x\t"<<"y\t"<<"Temp"<<" \  n"  ; 

double  outysize,  outxsize,  outzsize; 

outysize  = (outymax  -  outymin) /outynum; 

outxsize  = (outxmax  -  outxmin) /outxnum; 

outzsize  = (outzmax  -  outzmin) /outznum; 

//using  dummy  variables  in  case  one  wants  to  use  bigtherm2 (void)  afterwards 
for (int  i  =  0;  i<tnum;  i++) { 
t  =  tmin  +  i  *  tsize; 
for (int  ii  =  0;  ii<outznum;  ii++) { 

z  =  outzmin  +  ii  *  outzsize; 
for (int  iii  =  0;  iii<outxnum;  iii++) { 
x  =  outxmin  +  iii*outxsize; 
for (int  vi  =  0;  vi<outynum;  vi++)  { 

y  =  outymin  +  vi*outysize; 
double  thing  =  BigThermTempAt (x, y, z, t) ; 


writef ile<<t<<"\t"<<z<<"\t"<<x<<" \t"<<y<<"\t"<<thing<<" \n" ; 

} 

} 

writef ile<<" \n" ; 

} 

} 


} 

void  takatadamage (char*  conf igf ilename) { 


int  DamFlag; 

int  xnum,  ynum,  znum,  tnum; 

double  Pi,  BoxLength,  BoxDepth,  a,  MelDens; 
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double  xmin,  xmax,  xsize; 
double  ymin,  ymax,  ysize; 
double  zmin,  zmax,  zsize; 
double  tmin,  tmax,  tsize; 
vector<  double  >  Temp; 
double  ell,  c21,  cl2,  c22; 
double  irr , pulsedur , E ; 

double  x,  y,  z,  t,  Tbody,  rnum,  absorb; 
double  TempA; 

double  spotsize, spotsizex, spotsizey; 

float  Terml,Term2; 

float  Suml,  Sum2; 

float  Integrand; 

float  Damage; 

int  old; 

string  inputfile; 
string  outputfile; 

boost: : property_tree : :ptree  configfile; 

boost :: property  tree :: ini_parser :: read  ini (conf igf ilename,  configfile) 
string  OutputFile; 

// - initializing  the  output  file 

outputfile  =  conf igf ile . get<string> ( "Bigtherm. DamageOutputFile" ) ; 
f stream  out (outputfile . c_str ( ) ) ; 
out . open (outputfile . c_str ( ) )  ; 

out<<"t\t"<<"x\t"<<"y\t"<<" z\t"<<"Damage\t"<<endl ; 


//! - Initializing  constants: 

Tbody=37 . 0 ; 

DamFlag=0 ; 
cl 1=1 4  9 . 0 ; 
cl2=50000 . 0 ; 
c21=242 . 0; 
c22=80000 . 0 ; 

// - Reading  input  variables - 

BoxLength=conf igf ile . get<double> ( "Bigtherm. BoxLength" ) ; 
BoxDepth=conf igf ile . get<double> ( "Bigtherm. BoxDepth" ) ; 
xnum=conf igf ile . get<int> ( "Bigtherm. xnum" ) ; 
xmin=conf igf ile . get<double> ( "Bigtherm. xmin" ) ; 
xmax=conf igf ile . get<double> ( "Bigtherm. xmax" ) ; 
ynum=conf igf ile . get<int> ( "Bigtherm. ynum" ) ; 
ymin=conf igf ile . get<double> ( "Bigtherm. ymin" ) ; 
ymax=conf igf ile . get<double> ( "Bigtherm. ymax" ) ; 
znum=conf igf ile . get<int> ( "Bigtherm. znum" ) ; 
zmin=conf igf ile . get<double> ( "Bigtherm. zmin" ) ; 
zmax=conf igf ile . get<double> ( "Bigtherm. zmax" ) ; 
tnum=conf igf ile . get<int> ( "Bigtherm. tnum" ) ; 
tmin=conf igf ile . get<double> ( "Bigtherm. tmin" ) ; 
tmax=conf igf ile . get<double> ( "Bigtherm. tmax" ) ; 

MelDens=conf igf ile . get<double> ( "Bigtherm. MelDens" ) ; 
rnum=conf igf ile . get<int> ( "Bigtherm. rnum" ) ; 
old=conf igf ile . get<int> ( "Bigtherm. old" ) ; 
absorb=conf igf ile . get<double> ( "Bigtherm. absorb" ) ; 
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//! 


Calculating  necessary  parameters 


xsize  =  (xmax  -  xmin) /xnum; 

ysize  =  (ymax  -  ymin)/ynum; 

zsize  =  (zmax  -  zmin) /znum; 

tsize  =  (tmax  -  tmin)/tnum; 

/* 

I  'k-k'k-k'k'k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k-k-k'k-k'k'k'k'k-k'k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k-k'k'k'k'k-k'k-k'k-k'k'k'k-k'k'k'k'k'k'k'k-k'k-k 

!  Calculate  Damage (x, y, z, t) 

I  ■k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k'k'k'k-k-k-k'k'k-k-k'k'k-k'k'k'k-k-k'k'k-k-k'k'k-k'k'k 
*  / 

//  initializing  input  file 

inputfile  =  conf igf ile . get<string> ( "Bigtherm. tempf ile" ) ; 
if stream  damagef ile ( inputfile . c  str () ) ; 
double  dummy; 
string  line; 
int  dam; 
if (  old  ==0 ) 

getline (damagef ile,  line); 

//! - Stepping  in  x  direction - 

for(int  i=0 ; i<xnum+l ; i++) { 
x  =  xmin  +  i  *  xsize; 

//! - Stepping  in  y  direction - 

for (int  ii  =  0;  ii<ynum+l;  ii++) { 
y  =  ymin  +  ii  *  ysize; 

//! - Stepping  in  z  direction - 

for(int  iii  =  0;  iii<znum+l;  iii++) { 
z  =  zmin  +  iii  *  zsize; 

//! - Stepping  in  time - 

for (int  iv  =  0;  iv<tnum+l;  iv++) { 
t  =  tmin  +  iv  *  tsize; 

if (old  ==1 ) { 

damagefile  >>  dummy; 

Temp. push  back(dummy); 

} 

if (old  ==  0 ) { 

damagefile 

>>dummy>>dummy>>dummy>>dummy>>dummy; 

Temp. push  back(dummy); 

} 

} 

//! - Simpson's  Rule  for  Damage  Integral - 

Suml=0 . 0 ; 

Sum2=0 . 0 ; 
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if (Tbody  +  Temp[0]  <  50.0) 


Terml=exp (ell  -  cl2/ (273 

+ 

Tbody 

+ 

Temp 

!0] )  )  ; 

else 

Terml=exp (c21  -  c22/ (273 

+ 

Tbody 

+ 

Temp 

[0]  )  )  ; 

if (Tbody  +  Temp[tnum-1]  <  50.0 

) 

Term2=exp (ell  -  cl2/ (273 

+ 

Tbody 

+ 

Temp 

[tnum-1 ] 

else 

Term2=exp (c21  -  c22/ (273 

+ 

Tbody 

+ 

Temp 

[tnum-1 ] 

for (int  iv  =1;  iv<tnum;  iv+=2) 

{ 

if  ((Tbody  +  Temp[iv])  < 

50 

1.0) 

Integrand=exp (ell  -  (cl2/ (273  +  Tbody  +  Tempfiv]))); 

else 

Integrand=exp (c21  -  c22/ (273  +  Tbody  +  Temp[iv])); 
Suml=Suml+ Integrand; 

} 

Suml*=4 . 0 ; 

for(int  iv  =2;  iv<tnum-l;  iv  +=2) { 
if  ((Tbody  +  Temp[0])  <  50.0) 

Integrand  =  exp(cll  -  (cl2/ (273  +  Tbody  +  Temp [iv] ) ) ) ; 

else 

Integrand=exp (c21  -  c22/ (273  +  Tbody  +  Temp[iv])); 
Sum2+= Integrand; 

} 

Sum2*=2 . 0 ; 


Damage=tsize* (Terml  +  Term2  +  Suml  +  Sum2)/3.0; 

if (Damage  >=  1.0) 
dam  =  1 ; 

else 

dam  =  0 ; 


cout«x<<"\t"<<y«"\t"<<z«"\t"<<Damage<<endl; 


} 

} 

void 


} 


if  (Damage  >=1.0) 
DamFlag=l ; 
Temp . clear ( ) ; 


out . close ( ) ; 
damagef ile . close ( )  ; 

BigTherm: : TakataDamage (void) { 


int  DamFlag; 
vector<  double  >  Temp; 
double  ell,  c21,  cl2,  c22; 
double  x,  y,  z,  t; 
double  TempA; 

double  spotsize, spotsizex, spotsizey; 
float  Terml, Term2; 
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float  Suml,  Sum2; 
float  Integrand; 
float  Damage; 
int  old; 

// - initializing  the  output  file 

of stream  DamageOut (DamageOutputFile . c_str ( ) , ios : : out) ; 
of stream  TempOut (TempOutputFile . c_str ( ) , ios : : out) ; 

DamageOut  <<  "x\t"  <<  "y\t"  <<  "z\t"  <<"Damage\t"  «endl; 
if (Damagelnfo  ==  1){ 

TempOut<<"t\t"<<"x\t"<<"y\t"<<" z\t"<<"Temp\t"<<endl ; 

} 

//! - Initializing  constants: 

Tbody=37 . 0 ; 

DamFlag=0 ; 
cll=149 . 0; 
cl2=50000 . 0 ; 
c21=242 . 0 ; 
c22  =  80000 . 0 ; 

//! - Calculating  necessary  parameters - 


xsize  = 

(xmax  - 

xmin) 

/ xnum; 

ysize  = 

(ymax  - 

ymin) 

/ynum; 

zsize  = 

(zmax  - 

zmin) 

/ znum; 

tsize  = 

(tmax  - 

tmin) 

/ tnum; 

/* 

I  ■k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k-k'k-k'k-k'k'k-k-k'k-k'k-k'k'k-k-k'k'k-k-k'k'k'k-k-k'k'k-k-k'k'k-k-k'k'k-k'k'k'k-k-k'k'k-k'k'k'k-k'k'k 

!  Calculate  Damage (x, y, z, t) 

I  : k-k'k-k'k-k-k-k'k-k'k-k'k-k-k-k'k'k'k-k'k-k'k-k'k-k'k-k'k-k-k-k'k-k'k-k'k-k'k-k'k-k'k-k-k'k-k'k-k'k-k'k-k'k-k-k-k'k-k'k-k'k'k'k-k'k-k'k-k'k'k 
*/ 

//  initializing  input  file 
double  dummy; 
int  dam; 

//! - Stepping  in  x  direction - 

for (int  i=0 ; i<xnum+l ; i++) { 
x  =  xmin  +  i  *  xsize; 

//! - Stepping  in  y  direction - 

for (int  ii  =  0;  ii<ynum+l;  ii++) { 
y  =  ymin  +  ii  *  ysize; 

//! - Stepping  in  z  direction - 

for (int  iii  =  0;  iii<znum+l;  iii++) { 
z  =  zmin  +  iii  *  zsize; 

//! - Stepping  in  time - 
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for(int  iv  =  0;  iv<tnum+l;  iv++) { 
t  =  tmin  +  iv  *  tsize; 
dummy  =  BigThermTempAt (x, y, z, t) ; 
Temp. push  back (dummy) ; 
if (Damagelnfo  ==  1) 

TempOut<<t<<"\t"<<z<<"\t"<<x<<"\t"<<y<<"\t"<<dummy<<endl  ; 

} 

//! - Simpson's  Rule  for  Damage  Integral - 

Suml=0 . 0 ; 

Sum2=0 . 0 ; 

if (Tbody  +  Temp[0]  <  50.0) 

Terml=exp (ell  -  cl2/(273  +  Tbody  +  TempfO])); 

else 

Terml=exp (c21  -  c22/(273  +  Tbody  +  TempfO])); 

if (Tbody  +  Temp[tnum-1]  <  50.0) 

Term2=exp (ell  -  cl2/(273  +  Tbody  +  Temp [tnum-1] ) ) ; 

else 

Term2=exp (c21  -  c22/(273  +  Tbody  +  Temp [tnum-1] )) ; 

for(int  iv  =1;  iv<tnum;  iv+=2) { 

if  ((Tbody  +  Temp[iv])  <  50.0) 

Integrand=exp (ell  -  (cl2/ (273  +  Tbody  +  Tempfiv]))); 

else 

Integrand=exp (c21  -  c22/ (273  +  Tbody  +  Tempfiv])); 

Suml  =  Suml+  Integrand; 

} 

Suml*=4 . 0; 

for(int  iv  =2;  iv<tnum-l;  iv  +=2) { 
if  ((Tbody  +  Temp[0])  <  50.0) 

Integrand  =  exp(cll  -  (cl2/ (273  +  Tbody  +  Tempfiv]))); 

else 

Integrand  =  exp(c21  -  c22/ (273  +  Tbody  +  Tempfiv])); 
Sum2 += Integrand ; 

} 

Sum2*=2 . 0 ; 

Damage=tsize* (Terml  +  Term2  +  Suml  +  Sum2)/3.0; 
if (Damage  >=  1.0) 
dam  =  1 ; 

else 

dam  =  0 ; 

DamageOut<<x<<" \ t"<<y<<" \ t"<<z<<" \ t"<<dam<<endl ; 
if  (Damage  >=1.0) 

DamFlag=l ; 

Temp . clear ( ) ; 

} 

} 

} 

DamageOut . close ( )  ; 

TempOut . close ( ) ; 


} 
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7.4  Input  File 


;  Text  preceeded  by  or  are  comments  and  will  not  be  read  by 

;  boost_property_tree  parser.  This  example  will  generate  a  temperature 

;  profile  and  damage  spot  from  a  given  seed. 

[Bigtherm] 

BoxLength=  0.0040 

BoxDepth  =  0.0014 

GranularRadius=  0.0001 

xnum=40 

xmin=-0 .0020 

xmax=0 .0020 

ynum=40 

ymin=-0 .0020 

ymax=0 .0020 

znum=2 

zmin=-0 .0007 
zmax=0 .0007 
tnum=20 
tmin=0 . 0 
tmax=0 . 002 
Melnum= 

MelDens=8 . 33e9 
rnum=500 
absorb=4000 . 0 

; 0.0040  0.0014  .0001  #  BoxLength,  BoxDepth,  Granular  Radius (a) 

;40  -.0020  .0020  #  Number  of  x  steps (xnum) ,  xmin,  xmax 

;40  -.0020  .0020  #  Number  of  y  steps (ynum) ,  ymin,  ymax 

;2  -.0007  .0007  #  Number  of  z  steps (znum),  zmin,  zmax 

;20  0.0  .002  #  Number  of  time  steps (tnum) ,  tmin,  tmax 

;8.33D9  500  4000.0  #  MelDens,  rnum,  absorb 


Comments  here  will  not  be  read  by  C++ 

1)  Input  all  distances  in  cm  and  all  times  in  sec. 

Melanin  density (MelDens )  in  granules/cm**3 . 

Melanin  absorption  coefficient (absorb)  in  cm**-l. 

All  calculations  are  done  in  cgs  units.  Some  energies  are  input  in 

Joules;  but  these  are  converted  to  calories  before  calculation. 

2)  Melanin  Parameters 

A.  Melanosome  (or  granular)  radius  is  roughly  1  micron.  This  is 
a  variable  input  parameter. 

B.  Melanosomes  are  distributed  through  a  rectangular  volume  whose 
x,y  dimensions  are  BoxLength  and  z  dimension  is  BoxDepth. 
Granules  are  randomly  distributed  through  this  space  using  a 
Cartesian  coordinate  system  with  origin  at  the  geometric  center 
of  the  volume.  Granular  volumes  do  not  overlap. 

C.  The  volume  of  an  RPE  cell  is  independent  of  the  melanin  dist. 
volume.  It  is  estimated  to  be  20  x  20  microns  cross  section (x,y) 
and  15  microns  depth (z)  =  6  x  10**-9  cm**3.  The  portion  of  the 
RPE  cell  volume  that  contains  melanosomes  is  Vm  =  20  microns  x 
20  microns  x  BoxDepth. 
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cell . 


D.  The  melanin  number  density,  Nm,  is  the  number  of  granules/RPE 


using : 


gives 


The  melanin  volume  density,  MelDens,  may  be  computed  from  Nm 

MelDens (granules/cm**3)  =  Nm (granules/cell ) /Vm (cm**3/cell )  . 

For  example,  BoxDepth  =  15  microns  and  Nm  =  120  granules/cell 

Vm  =  6E-9  cm**3/cell  and  MelDens  =  2E10  granules/cm**3 . 

E.  The  number  of  granules  in  the  whole  dist.  vol. (melnum)  is: 
melnum  =  MelDens  *  BoxLength**2  *  BoxDepth 
To  avoid  excessive  cpu  time,  this  number  should  not  exceed  20,000. 
(The  code  is  currently  hard-wired  for  arrays  of  dimension  20,000) . 


3)  The  single  granule  temperature  distribution,  T(r,t)  is  calculated 
using  a  spherical  coordinate  system  independent  of  the  Cartesian 
system  used  for  melanin  dist.  The  calculation  is  radially  symmetric 
and  is  done  for  rnum  radial  steps  and  tnum  time  steps.  Note  that 
T(r,t)  is  analytically  calculated  and  is  not  the  product  of  a  temporal 
summation  or  a  spatial  stepping  algorithm.  Each  value  can  thus  be 
calculated  independently  of  values  at  other  positions  and  times. 


;  4)  The  final  temperature  distribution  in  the  retina  is  calculated  using 

;  the  same  Cartesian  coordinate  system  as  for  melanin  dist.  The  value 

;  T(x,y,z,t)  at  a  certain  arbitrary  position  and  time  is  the  linear 

;  superposition  of  the  contributions  of  all  the  granules  T (rl , t) +T (r2 , t) 

;  +T(r3,t)+...  at  that  time,  where  rl , r2 , r3 , etc .  are  the  radial 

distances 

;  from  each  granule  to  (x,y,z).  The  value  T(x,y,z,t),  like  T(r,t),  can 

;  be  calculated  independently  for  each  position  and  time. 


;  5)  Since  each  T(x,y,z,t)  is  independently  calculated,  the  space  and  time 

;  grid  chosen  for  the  temperature  calculation  is  in  theory  arbitrary. 

;  Note  that  the  number  of  steps  in  each  dimension  is 

xnum, ynum, znum, tnum; 

;  whereas  the  actual  number  of  grid  points  is  xnum+1 , ynum+1 , znum+1 ,  and 

;  tnum+1 .  The  do  loops  in  the  code  go  from  0  to  xnum, ynum, etc .  and  not 

;  from  1  to  xnum+1, etc. 

;  6)  The  number  of  x,y,&  z  grid  points  is  only  constrained  by  the  spatial 

;  resolution  needed  for  visualization.  Max  and  min  values  may  be  inside 

;  or  outside  the  melanin  volume.  They  need  not  be  symmetric  about  the 

;  origin.  They  are  not  tied  to  the  maximum  value  of  r  used  in 

calculating  ;  T(r,t)  for  one  granule.  For  r>rmax  the  contribution  of 

that  granule  is  ;  simply  0. 

r 

;  7)  Typically  xnum  &  ynum  should  be  set  to  give  at  least  1  step/micron.  A 

;  minimum  of  znum+1  =  3  should  be  used  where  the  3  z  slices  represent 

the 

;  top,  middle,  and  bottom  layers  of  the  RPE .  More  z  steps  may  be  needed 

;  to  accurately  construct  a  side  view  of  the  damage  volume  or  to 

;  accurately  account  for  true  shadowing. 

r 

;  8)  The  number  of  time  steps  tnum  is  constrained  by  the  requirements  of 

;  the  damage  calculation.  Temporal  resolution  must  be  great  enough  to 

;  give  convergance  of  the  Arrenhius  integral  over  the  time-temperature 

;  history.  Typically  tnum  of  about  20-40  is  sufficient.  Unlike  the 

;  spatial  grid  the  temporal  calculation  is  limited  to  those  time  values 

;  for  which  T(r,t)  has  previously  been  calculated. 
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;  9)  The  minmum,  tmin,  is  always  0  (before  the  pulse) .  The  max  value  can  be 

;  inside  or  outside  the  pulse  interaction  time,  if  the  only  thing 

desired 

;  is  the  temperature  distribution.  For  damage  calculations  to  be 

accurate 

;  however,  tmax  must  be  large  enough  that  the  tissue  has  cooled  to 

normal 

;  (37  C)  .  For  most  cases  the  pulsewidth  plus  2  msec  should  be 

sufficient . 

;  Very  long  pulses  (on  the  order  of  seconds)  or  very  large  spots  may 

;  reqiuire  more  than  2  msec  cooling  time  and  adjustments  may  be  needed. 

;  If  tmax  is  large,  the  energy  conservation  report  will  likely  show 

;  that  much  energy  has  been  lost.  This  simply  means  that  energy  has 

;  had  time  to  move  out  of  the  RPE  and  into  the  surrounding  tissue. 

r 

;  10)  It  is  not  a  good  idea  to  run  this  model  for  times,  t  <  1  microsec. 

;  (The  pulse  duration  can  be  as  short  as  you  desire,  but  in  the 

;  ultrashort  case,  don't  ask  for  a  temperature  profile  until  well  after 

;  the  pulse  has  expired.  We  do  not  expect  that  thermal  mechanisms  are 

;  relevant  for  ultrashort  anyway.) 

f 

;  11) The  Arrhenius  Integral  damage  calculations  are  done  using  a  Simpson's 

;  Rule  Integrater.  The  coefficients  are  from  either  Takata  or  Welch. 

r 

;  12) Image  diameter-  The  BoxLength  should  be  greater  than  or  equal  to  the 

;  transverse  beam  diameter  of  the  laser;  i.e.,  the  beam  should  encounter 

;  melanin  everywhere  it  falls  on  the  retina.  It  is  probably  best  if 

twice 

;  the  laser  diameter  is  used,  so  that  all  the  beam  is  covered  out  to  low 

;  intensities;  however,  this  not  strictly  necessary  for  validity. 

f 

;  13) Run  Time  -  There  are  basically  four  parts  to  the  thermal  code: 

r 

;  i.  Random  granule  placement  -  varies  with  melnum. 

;  ii.  Depth  averaging  or  shadowing  -  varies  with  melnum. 

;  iii.  T(r,t)  single  granule  calc.  -  varies  with  tnum. 

;  iv.  T(x,y,z,t)  mult.  gran,  calc.-  varies  with  xnum,  ynum,  znum,  tnum, 

;  melnum 

r 

;  Parts  i,ii,  and  iii  do  not  take  a  high  percentage  of  cpu  time  even  at 

;  fairly  high  values  of  melnum  and  tnum.  Almost  all  cpu  time  is  spent  on 

;  T(x,y, z,t)  which  is  linear  in  xnum, ynum, znum, tnum,  and  melnum. 


; inf ormation  formerly  contained  in  the  inp2  files 
irr=3 . Oe-5 

; (irr) Corneal  Fluence  in  J/cm**2. 
pulsedur=3 . Oe-6 

;  (pulsedur) Pulse  Duration  in  sec 
alpha=0. 00132798 

;  (alpha) Thermal  Diffusivity  of  Medium  in  cm**2/sec. 
cond=0 .00556 

;  (cond) Thermal  Conductivity  of  Medium  in  J/cm-C-sec. 
corspot=0 . 7 

;  (corspot) Beam  Diameter  at  cornea  in  cm. 
spotsize=0 .0022 
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;  (spotsize) Retinal  Image  Diameter  at  l/e**2  points  in  cm. 
trans=l . 0 

;  (trans) Fraction  of  pulse  energy  transmitted  to  retina, 
improf =0 

;  ( improf ) Image  Type [0=Gaussian,  3=Top  Hat,  4=Annular] . 

SeedFlag=l 

;  SeedFlag  [0=Fixed  Seed,  1  =  Random  Seed  from  time()  function] . 

;  1)  Input  all  distances  in  cm  and  all  times  in  sec. 

;  All  calculations  are  done  in  cgs  units.  Some  energies  are  input  in 

;  Joules;  but  these  are  converted  to  calories  before  calculation. 

r 

;  2)  Corneal  beam  diameter  is  used  to  compute  the  Total  Interocular 

;  Energy  (multiplying  area  by  Corneal  Fluence) .  Therefore,  if 

;  the  beam  is  larger  than  the  pupil,  the  pupil  diameter  should  be  used. 

r 

;  3)  The  average  focusing  power  of  the  eye  can  be  determined  by: 

;  Retinal  Fluence  /  Corneal  Fluence  =  (Corneal  Diam  /  Retinal  Diam)**2. 

;  For  10**5,  you  might  choose  Corneal  Diam  =  .7  and  Retinal  Diam  = 

.0022  . 

r 

;  4)  Fraction  of  pulse  energy  transmitted  is  a  function  of  wavelength  and 

;  may  be  obtained  from  the  literature.  (See  Maher) 

r 

;  6)  Values  of  thermal  diffusivity  and  conductivity  for  water: 

;  0.00132798  cm**2/sec  and  0.00556  J/cm-C-sec 

;file  for  explicit  placement  of  melinen  granules 
MelPlacementFile=ExampleMelPlacementFile . txt 

TempOutputFile=ExampleOutputTempFile . temp 

; name  for  output  temperature  file,  the  file  extention  is  arbitrary 

;  output  will  be  formated  time,  x,  y,  z,  temperature,  units  are  the  same  as 

above 

DamageOutputFile=ExampleOutputDamageFile . dam 

; name  for  output  damage  file,  the  file  extention  is  arbitrary, 

; output  will  be  x,  y,  z  then  either  a  0  or  1 

;0  indicates  no  thermal  damage  has  been  done,  1  indicates  damage  has  been 
done 

PulseNum=l 

; The  number  of  pulses  long  an  exposure  is,  defauts  to  1  if  none  is  given 
PulseSep=0 .001 

; the  time  in  seconds  between  each  pulse, 
tempf ile=ExampleInputTempFile . temp 

;If  you're  using  the  simple  re-implementation  of  the  takata  model,  this  is 

the  temperature  file  you  need  for  that. 

old=0 

;this  was  an  option  that's  only  needed  in  the  re-implementation  of  the  takata 
model,  then  you  need  to  set  this  flag  to  1  in  order  to  use  the  old  fortran 
temp  files. 

Damage Inf o=l 

;  set  this  flag  to  1  if  you  want  the  temperature  file  generated  for  the  damage 
to  be  output  as  well 
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